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Abstract 

We analyze a class of network motifs in which a short, two- node positive feed- 
back motif is inserted in a three-node negative feedback loop. We demonstrate 
that such networks can undergo a bifurcation to a state where a stable fixed 
point and a stable limit cycle coexist. At the bifurcation point the period of 
the oscillations diverges. Further, intrinsic noise can make the system switch 
between oscillatory state and the stationary state spontaneously. We find that 
this switching also occurs in previous models of circadian clocks that use this 
combination of positive and negative feedback. Our results suggest that real- 
life circadian systems may need specific regulation to prevent or minimize such 
switching events. 

Keywords: negative feedback, positive feedback, oscillation, noise, circadian 
rhythm 



1. Introduction 

Genetic regulatory networks exhibit a wide range of oscillatory phenomena, 
ranging from the fast oscillations in calcium 1J to the 24 hour cycle of circadian 
clocks [21 [3] . The occurrence of oscillations is generally caused by the presence 
of a negative feedback loop in the regulatory network [?J [5] . From a theoretical 
point of view, negative feedback can cause a Hopf bifurcation and thus generate 
a transition between a stable fixed point, corresponding to homeostasis, and an 
attracting limit cycle, corresponding to oscillations. 
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However, in real regulatory networks, the loop causing oscillations is usually 
embedded in a larger network including multiple positive and negative feed- 
backs. In some cases, these additional loops have a demonstrable biological 
function, for example in giving tunability to the oscillation period [6 or in 
stabilizing the period of circadian clocks in the presence of temperature fluctua- 
tions or molecular noise [7j. In general, one expects that multiple feedback loops 
could lead to non-trivial behaviors from the viewpoint of bifurcation theory. The 
dynamics could become even richer when the noise induced by stochastic gene 
expression is taken into account. As we discuss later, circadian clocks are typical 
examples of genetic circuits where it may be important to understanding these 
effects. 

In this paper, we present and analyze a class of network motifs in which a 
two-node positive feedback motif is inserted into a three-node negative feedback 
loop, as represented in Fig. [T] In section 2.1, we show that a deterministic dy- 
namical models of the simplest of such networks, in a suitable parameter range, 
exhibits co-existence of a stable fixed point and a stable limit cycle. We explain 
this behavior in terms of a saddle- node separatrix-loop bifurcation [5] , and show 
that it results in a diverging oscillation period close to the bifurcation point. 
In section 2.2, we use stochastic simulations using the Gillespie algorithm [9] 
to demonstrate that the noise can make the system switch between oscillatory 
state and the stationary state. We show that similar behaviour occurs in more 
realistic models of circadian clocks that contain the same combination of pos- 
itive and negative feedback loops. Section 3 summarizes our thoughts on the 
relevance of these results for the behaviour of circadian clocks. 

2. Results 

2.1. Dynamics of the Network Motif 

We study the class of genetic networks represented in Fig. [I] In each of 
the four networks, a positive feedback between node 1 and node 2 can give rise 
to a bi-stable switch. When node 3 is introduced, node 1, 2, 3 together form 
a negative feedback loop. The negative feedback loop tends to destabilize one 
of the stable fixed points of the switch. The simplest motifs that exhibit such 
"frustrated bistability" are studied in ref. [10]. Here, we study slightly larger 
networks which allow for more intricate dynamics, in particular a scenario where 
a stable limit cycle emerges around the unstable fixed point while the other 
stable fixed point remains unchanged. 

We shall first focus our discussion on the network (a) in Fig [I] for which we 
write the following dynamical equations: 
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Here, £1,2,3 are the concentrations of the proteins associated with the three 
nodes, a is the strength of the three inhibitory regulations, h is the Hill coeffi- 
cient, c is a constant source term for each node, and 7 is the degradation rate 
for each protein (we assume that all three proteins are stable, and therefore 
their degradation rate is determined by the cell division time). j3 is the con- 
trol parameter to adjust the inhibition from node 2 to node 1. To simplify the 
equations, we introduce dimensionless parameters and variables X\ 23 = X1 L 2 ' 8 , 
T = 7 t,A=%,B=^: 

^ = B + A 1 h hm: -* W 

dr 1 + X 3 h l + {^) h 

f - B + ^W"^ <6 » 

We study this system of equations for parameter values of h = 3,B = 
0.1, A — 5, and vary (3 as a control parameter, which changes the strength of the 
positive feedback relative to the negative feedback. We found three bifurcation 
points, at (3 = A 3 ~ 2.03, /3 = A 2 « 1.923, and $ = Ai « 1.45. When /3 > A 3 
(Fig[2j\_), only one unstable fixed point and one stable limit cycle are found in 
the phase space. Here, the stable limit cycle is a global attractor. At j3 = A3, 
a saddle-node bifurcation occurs, where a stable fixed point and a saddle node 
emerge. 

When A2 < (3 < A3 (Figj2j3), we find one stable limit cycle and three fixed 
points - one stable and two unstable. Within this region, the stable limit cycle 
is not the global attractor. Depending on the initial condition, the system may 
either reach the limit cycle or the stable fixed point. In the parameter range 
we studied, the volumes of the two basins of attraction are both non-negligible, 
separated by the surface shown in Fig [3] 

Approaching the critical point A2 , the stable limit cycle approaches the sad- 
dle point, and eventually at (3 = A 2 (Fig{2]C), they generate a homoclinic cycle. 
When (3 is near A2, the period of oscillation increases dramatically, and even- 
tually diverges at the critical point f3 = \ 2 (Fig [4} as typical for homoclinic 
cycles. Such a bifurcation is classified as saddle-node separatrix-loop bifurca- 
tions [SI [11] and is a robust bifurcation scenario in a phase space of dimension 
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3 or more |5]. It corresponds to an attractor crisis, so that when (3 is reduced 
further, the stable limit cycle disappears. When Ai < ft < A 2 (Fig{2p), three 
fixed points still exist in the phase space: a stable fixed point, a saddle node and 
an unstable fixed point with complex eigenvalues. The stable fixed point is the 
global attractor. At /? = Ai, the saddle node and the unstable fixed point collide 
and disappear after a saddle-node bifurcation. So, when j3 < Ai (Figj^), there 
is only one fixed point, which is stable and is a global attractor. 

The location and nature of the fixed points can be better understood by 
a graphical study of the intersections of the nullclines. We set d/dr terms 
in equations |4|-([6]) to zero, and rearrange the resulting algebraic relations to 
express Xi and X3 in terms of X\ . Then, using this to eliminate X2 and X3 in 
equation Q yields: 

= TT7 ^ igHi =3 - + b-x 1 . (7) 

Fig [2] right column shows the right hand side of Eq. Q as the parameter 
j3 is varied. When j3 > A3 (Fig(2j^) , the function has one zero, corresponding 
to a unique unstable fixed point in phase space. When A 2 < /3 < A3 (Fig[2j3), 
/3 = A 2 (Figj2p) and Ai < j3 < A 2 (Fig{2p), the function has two additional 
zeroes, indicating two more fixed points. From the eigen values of all three 
fixed points, one can show that only one of them (the circular one) is stable. 
Finally, when /? < Ai (Figj^H), the function has one zero again. Two fixed 
points disappear and only one stable fixed point remains. 

We checked numerically the parameter range where the same bifurcations 
and qualitatively the same phase space portrait can be obtained. For A = 5 and 
B = 0.1, h>,3 is required to see the same behavior. For h — 3 and B = 0.1, 
A > 3.5 is required. The behavior is found to be insensitive to the value of B: 
for h = 3 and A = 5, the behavior was unchanged for < B < 100. The same 
bifurcations can also be found in all the other 3 motifs listed in Fig [l] Namely, 
the observed sequence of bifurcations is a robust feature of such motifs. 

2.2. Stochastic Dynamics 

In this section, we investigate the effect of the intrinsic noise due to the 
discrete nature of molecular reactions in such motifs. We use the Gillespie 
algorithm [5] for stochastic simulations of the dynamical system specified by 
equations ([l])-([3]). We denote the copy number of molecules of species i as iVj. 
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The allowed transitions, along with their kinetic rates, are: 
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To control the noise, we change the volume of the system V, which changes 
the production rates of 2Vj, but leaves the average concentration = N/V 
constant, as long as the values of k, 7, h, c, a, and /? are unchanged. The larger 
the value of V, and therefore the larger the copy numbers Ni, the smaller the 
noise (the relative fluctuations in Ni). Note that we do not explicitly consider 
processes like mRNA production, binding of transcription factors, etc. Inclusion 
of these steps can increase the noise in the system further [12] . 

Figure [5] shows the concentration N/V vs. time for a stochastic simulation 
in dimensionless units with h = 3, B = 0.1, A = 5, j3 = 2. We convert numbers 
so that the dimensionless concentration Xi = 1 corresponds to one molecule 
when volume V = 1, and simulated (A) V = 1000 and (B) V = 100. We can 
clearly see switching between the oscillatory state (the stable limit cycle) and 
and the steady state (the stable fixed point), and this switching happens more 
often for smaller V, i.e., for larger noise. 

To quantify this switching behavior, we measured the average switching 
time from the oscillatory state to the steady state and vice versa as a function 
of the system volume V (Fig{6]). Because of the noisy dynamics, switching is 
determined by using two thresholds for the distance from the stable fixed point, 
Si and S2 > Si : We define a switching event from the steady regime to osillatory 
regime when when the distance exceeds 5*2 , while the reverse switching happens 
when the distance becomes smaller than Si. Therefore, when the distance is 
between Si and 5*2, there is a history dependence in which regime the state 
belongs to. For this parameter set, the oscillation period is ~8.15, thus we can 
see that the system with hundreds of molecules (corresponds to V w 100) can 
still cause frequent switching, on the order of once in every 10 oscillations. The 
switching rate decreases with V as expected. For large enough V, the switching 
rate decreases exponentially with V. 

Finally, in order to see whether such switching behavior can be relevant for 
real biological systems, we study effect of noise on the more realistic model 
for Drosophila circadian rhythms described in ref. [131 [T3] . The deterministic 
version of this model has a combination of several positive and negative feedback 
loops and exhibits the coexistence of a stable fixed point and a stable limit cycle 
[13] . We simulated the stochastic version of this model with parameters used in 
ref. [14j . A detailed description of the model and parameters are given in the 
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appendix. We observe again the switching between the oscillatory state and the 
steady state due to the noise (Fig [7| . The switching is quite often compared to 
the oscillation frequency for the noise level expected for an average cell volume, 
i.e. 10 3 ^m (Fig. Up), where 1 [nM] corresponds to about 600 moleculeaM 




3. Summary and Discussion 

We have analyzed a class of network motifs consisting of a two-node positive 
feedback inserted into a three-node negative feedback loop. We demonstrated 
that a stable fixed point and a stable limit cycle can co-exist in this class of mo- 
tifs. As parameters are changed, the system undergoes a saddle- node separatrix- 
loop bifurcation, with a diverging oscillation period as the system approaches 
the bifurcation. The location and nature of the fixed points were investigated in 
detail by outlining the intersections of the nullclines of the three variables. We 
then studied the effect of intrinsic noise to the motif, in the parameter regime 
where a fixed point and stable limit cycle co-exist. We showed that stochastic 
switching between the two happens with a rate that decreases with decreas- 
ing the level of the noise. Actually, our results complement the study in [25] . 
where similar motifs for biochemical systems are studied focusing on the how a 
negative feedback perturbs the switch by positive feedback loops. 

We further showed that similar behaviour is observed in a more complex 
model of the Drosophila circadian clock. This switching behaviour was not 
reported in a study of the effect of noise in a simplified version of the detailed 
model [7] . The studies on the Drosophila model [121 HH [7] focused instead on 
the singular behaviour of the real circadian clock, where a single pulse of light 
can cause long-term suppression of the circadian rhythm [21 [15j [13] , The models 
explain this as a switch from the limit cycle to the stable fixed point, caused 
by a short external perturbation (the pulse of light) to some parameter values. 
Recent research suggests an alternate possibility, where the singular behaviour 
is caused by desynchronization of the clocks rather than stopping individual 
oscillations [TB] , 

Our analysis of simple motif combining positive and negative feedback demon- 
strates that intrinsic noise, and not just external perturbations, can also cause 
switching where a limit cycle and a stable fixed point coexist in the phase space. 
As such repeated switching would disrupt the circadian rhythm, we can predict 
that specific regulatory mechanism must exist in real circadian clocks to sup- 
press it. It would be interesting to explore the space of small network motifs 
to understand what mechanisms could implement this kind of suppression most 
effectively. 



Calculated based on 1 nM ~ 6 X 10 23 X 10 -9 molecules per liter. 
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Appendix: The model for Drosophila circadian rhythm 

We studied the stochastic version of the model for Drosophila circadian 
rhythm given in ref [21 [17]. The summary of the reactions in the model with de- 
terministic equations of the model are shown in Fig. [8] with parameters given in 
the caption. We converted these equations into stochastic form using the Gille- 
spie algorithm 9J, as has been done to convert the simple motifs eqs. 0-([3]) to 
the stochastic version (|8])-(13). The concentrations are converted to the number 



of molecules based on the typical cell volume of drosophila, about lOOO^im 3 , 
which means 1 nM corresponds to about 600 molecules. Figure [7p is based on 
this conversion. For the case where the cell volume is 10 (100) fold bigger, the 
copy number is also converted to be 10 (100) fold bigger, which corresponds to 
the simulations shown in Fig. (7^5 (A). 
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Figure Captions 



Figure 1: Four patterns of motifs, each consisting a switch and a negative feedback loop. Node 
1 and node 2 form a switch. Node 1,2,3 form a negative feedback loop. 



Figure 2: The left column is the real 3 dimensional phase space. The middle column is 2 
dimensional illustration of the real phase space. The right most column shows the right hand 
side of Eq. |7|. Parameters are set to be h = 3, B = 0.1, and A = 5. (A) A3 < f3, with one 
stable limit cycle (line) and one unstable fixed point, (plot with f3 = 2.2) (B) A2 < /3 < A3, 
with three fixed points (an unstable fixed point with complex eigenvalue (square), saddle 
node (triangle) and a stable fixed point (circle)) and a stable limit cycle, (plot with j3 = 2) 
(C) fi = A2, the saddle node hits the limit cycle and forms a homoclinic orbit, (plot with 
fi = 1.923) (D) Ai < /3 < A 2 , there are only three fixed points, (plot with = 1.2) (E) fi < Ai, 
there is only one stable fixed point, (plot with (3 = 2) 



Figure 3: The stable limit cycle and the stable fixed point are shown by solid line and filled 
circle, respectively, for h = 3, B = 0.1, A = 5, and /3 = 2. The surface shows the boundary 
between the basin of attraction of the stable limit cycle and the basin of attraction of the 
stable fixed point. 
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Figure 4: Period of oscillations as a function of 0, with h = 3, B = 0.1, A = 5. The period 
diverges at j3 = A2 ~ 1.923 (indicated by a vertical line), where a homoclinic cycle is present. 
As (3 — ¥ 00 the repressor link from node 2 to node 1 vanishes. What is left is a standard 
repressilator and the period one obtains in that limit (equal to 3.7 for these parameter values, 
shown by a horizontal line) is the period of this respressilator. 



Figure 5: Time evolution of the concentrations for the stochastic simulation of the system 
with h = 3, B = 0.1, A = 5,/3 = 2. All the units are dimensionless, and concentrations are 
converted to the numbers so that Xi = 1 corresponds to one molecule when V = 1. (A) 
V = 1000 and (B) V = 100. 



Figure 6: Switching rates from the oscillatory state to the steady state (open circles) and 
from the steady state to the oscillatory state (filled circles), as a function of system volume 
V , for the same parameter values as in Fig. [5] The threshold values are set to be Si = 0.316, 
S 2 = 3.317. 



Figure 7: Intrinsic noise is introduced to the original model of Drosophila (131 114) . Similar 
phenomena are observed. When cell size is very large, namely about 100 times of a typical 
cell volume, i.e. 10 5 ^tm (A), the system would switch from oscillation state to steady state at 
a random moment due to noise. As the cell size goes down (10 4 /im in B, 10 3 /im in C), the 
switching becomes more and more frequent. A typical cell volumes for Drosophila is about 
10 3 /im (C), which exhibit very noisy dynamics. Detailed description of the model with the 
parameter set is given in appendix, Fig. [8] 



Figure 8: The model for Drosophila circadian rhythm given in ref |14II17| . Negative feedbacks 
are shown with dotted line, while positive feedback loops are shown with dashed line. The 
deterministic equation of the model is shown in the right panel, where variables are concen- 
trations of per (M p ) and tim (My) mRNAs, the PER and TIM with three phosphorylation 
levels of Po (T ), Pi (Ti), and P 2 (T 2 ), respectively, the PER-TIM complex C, and the nuclear 
form of the PER-TIM complex (Cm)- Parameters used are: n = 4, v sP = 1.1 nM/h, v a x = 
1 nM/h, v mP =1.0 nM/h, v mT = 0.7nM/h, v dP = 2.2nM/h, k sP = k sT = 0.9 /h, fei = 0.8 
/h, k 2 = 0.2 /h, k 3 = 1.2/(nM- h), fc 4 = 0.6 /h, K mP = K mT = 0.2 nM, K IP = K IT =lnM, 
K dP = A" [JT =0.2nM, K 1P = K 1T = K 2P = K 2 t = K ZP = K ZT = K iP = A' 4T =2nM, 
V 1P = Vi T =8nM/h, V 2P = V 2 T=lnM/h, V 3P = V 3T =8nM/h, V iP = V 4T = 1 nM/h, 
k d = k dC = k dN = 0.01 /h, V dT = 1.3nM/h. 
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Figure 3 
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Figure 5 
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Figure 6 
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Figure 7 
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Figure 8 
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